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Abstract 

Biological protein interactions networks such as signal transduction or gene transcription networks are often treated as 
modular, allowing motifs to be analyzed in isolation from the rest of the network. IVlodularity is also a key assumption in 
synthetic biology, where it is similarly expected that when network motifs are combined together, they do not lose their 
essential characteristics. However, the interactions that a network module has with downstream elements change the 
dynamical equations describing the upstream module and thus may change the dynamic and static properties of the 
upstream circuit even without explicit feedback. In this work we analyze the behavior of a ubiquitous motif in gene 
transcription and signal transduction circuits: the switch. We show that adding an additional downstream component to the 
simple genetic toggle switch changes its dynamical properties by changing the underlying potential energy landscape, and 
skewing it in favor of the unloaded side, and in some situations adding loads to the genetic switch can also abrogate 
bistable behavior. We find that an additional positive feedback motif found in naturally occurring toggle switches could 
tune the potential energy landscape in a desirable manner. We also analyze autocatalytic signal transduction switches and 
show that a ubiquitous positive feedback switch can lose its switch-like properties when connected to a downstream load. 
Our analysis underscores the necessity of incorporating the effects of downstream components when understanding the 
physics of biochemical network motifs, and raises the question as to how these effects are managed in real biological 
systems. This analysis is particularly important when scaling synthetic networks to more complex organisms. 
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Introduction 

A longstanding question about signal transduction and gene 
transcription networks is how modular are they. Here modularity 
means relative insulation of small subgraphs or motifs of the main 
network from each other [1]. This question is especially relevant 
for synthetic biology that aims to build artificial circuits from the 
bottom up [2] . It is also relevant for molecular biologists that aim 
to arrive at a quantitative understanding of a cellular decision, by, 
for example, isolating a crucial network module [3] . 

For synthetic biologists the challenge is now to move from 
simple network motifs such as pulse generators [4], genetic 
switches [5-8], logic gates [9,10], and oscillators [11-13] to more 
complicated networks combining multiple motifs and networks in 
more complex organisms. Novel applications currently being 
explored include plant biosensors [14], hazardous waste remedi- 
ation [15], clean fuel technology [16], and numerous medical 
applications [17-20]. Synthetic biologists hope to utilize biological 
modules in a manner similar to electrical circuit board compo- 
nents - plugging them together to attain a specific, and novel, 
function [21]. At the core of the concept of either breaking down 
complex biological systems into small modules, or even building 
complex systems from modules, is the belief that these modules will 



behave predictably in isolation and in connection. Recent 
theoretical and experimental work however [22-25] suggests that 
the functioning of modules may not be independent of the 
downstream components that they are connected to. Adding an 
additional binding reaction to the output of a gene regulatory 
network (or loading the network) may decrease system bandwidth 
[24] and substrate sequestration in covalent modification cycles 
may result in signaling delay [26]. In vitro studies find that there is 
significant load-induced modulation of the upstream module in an 
enzymatic signal transduction cascades [24]. Theoretical analysis 
has also shown that a load can change the fundamental properties 
of an oscillating circuit [27]. Thus understanding the effects of 
adding a load to the output of these technologically important 
network modules is required for a thorough understanding of the 
challenges of scaling up synthetic networks to higher levels of 
complexity. 

Loads could also have noteworthy unrecognized effects in 
natural systems. In fact aU natural systems have loads in some ways 
or the other. Motifs in signal transduction networks are connected 
directly to a transcriptional response, or to downstream proteins 
that may function as transcription factors or go on to activate tran- 
scription factors. Motifs in gene transcription networks have tran- 
scriptional outputs with protein domains that bind nonspecifically 
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Author Summary 

Cells rely on complex networks of protein-protein interac- 
tions in order to carry out life functions. Scientists believe 
that these networks are organized in a modular fashion; 
that is they are made up of functionally distinct parts like 
an electronic circuit. Modularity implies that just as we put 
together electronic parts to make an amplifier that we can 
use in many different circuits, we can put together 
biochemical reactions to make an amplifier, or a switch 
or an oscillator, which perform the same function in 
different organisms. This assumption is important in 
synthetic biology, where we engineer and assemble 
synthetic genetic circuits in living organisms in a modular 
fashion. We show that for important modules like genetic 
and signaling switches, the assumption of modularity has a 
crucial limitation. We show that if one simply connects a 
biological switch to another downstream circuit, the 
presence of the connection changes the operation of the 
switch, which in some cases may stop behaving like a 
switch. Our work underscores the importance of taking 
into account these downstream connections and suggests 
that natural systems may be balancing some of these 
components in order to ensure that despite diversity, 
modules continue to work in different systems with 
fidelity. 

and specifically to binding sites on the DNA, apart from inter- 
acting with other transcription factors. 

Circuits that function as switches play an important role in all 
biological signaling and gene transcription networks because they 
encode decisions. This change of state can be brought about by an 
external signal, or an internal accumulation of a protein, which 
can drive the system to a different steady state. Examples are the 
regulatory circuits for the cell cycle in yeast [28], mitogen-acti- 
vated protein kinase cascades in animal cells [29-3 1] , and the lysis- 
lysogeny switch in the X phage [32]. Since many small circuits can 
show this kind of behavior, switches are among the earliest and 
most well studied of protein interaction circuits [33]. The genetic 
toggle switch, which was one of the first two synthetic circuits 
constructed, is a well-known synthetic example [5]. Given the 
ubiquity and importance of switch-like motifs, it is important to 
understand how their function could be affected by binding down- 
stream partners. 

These reasons prompted our theoretical study of the behavior of 
a simple genetic toggle switch [5], a toggle switch with positive 
feedback as well as a common positive-feedback based switch 
involving Ras activation in lymphocytes [29,30] under a load on 
either one or both of its outputs. These circuits are shown in Fig. 1 
and described below. The simple toggle switch is a widely studied 
and emulated .synthetic network motif based on the mutual 
repression of two repressor proteins. However, naturally occurring 
toggle switches are often found connected to an additional positive 
autoregulatory component. For example in the competence system 
in B. subtilis, ComK represses the production of Rok and Rok 
represses the production of ComK; however ComK also has a 
strong positive feedback upon its own production [34]. Another 
example is found in the apoptosis network of many multicellular 
organisms, including mammals. Within the pathway controlling 
intrinsic apoptosis is a set of genes with double-negative repression, 
Casp3 and XIAP, again accompanied by positive autoregulation 
of Casp3 [35]. 

The Ras protein is a G-protein found on mammalian cellular 
membranes that is important in many cellular processes and is 
an upstream activator of the MAPK pathway. Ras goes from a 



GDP-bound inactive form to a GTP-bound active form, often in a 
digital manner [30], and previous studies in lymphocytes have 
shown that RasGDP is activated to RasGTP via a bistable switch 
that arises from a positive feedback loop on its own activation via 
SOS (Son of Sevenless) [30]. However the Ras switch very 
naturally has an associated load, since to transduce the cellular 
signals down along the MAPK/ERK pathway, RasGTP naturally 
binds to Raf kinase. Thus the Ras switch system contains all the 
elements we need to study the effects of adding loads to a bistable 
switch which is based on a positive feedback loop. 

Methods 

Genetic toggle switch 

The basic genetic toggle switch consists of two mutually repress- 
ing genes as shown in Fig. 1 along with an additional system to 
toggle the states. As shown in previous studies, with the right 
combination of parameters, the toggle switch will stay in one of 
two stable states, each characterized by a high concentration of 
one of the repressor proteins, and strong repression of the other. 
The toggle switch can now be induced to switch states using two 
possible strategies for inducing a transition: decrease the level of 
highly expressed protein [5,36] or increase the expression of one of 
the repressed proteins (Fig. 1) using an additional inducible system 
[36]. For a model which utilizes the latter protocol we obtain a 
system of four differential equations [36] after including a load. 
The load may be a protein, a small molecule or a binding site on 
DNA such that the bound complex prevents the repressor from 
binding to and repressing its conjugate promoter. In order to make 
the simplest and the most general model, we have assumed here 
that the repressors reversibly bind the load only in one copy. We 
assume that the total load Lit is a constant, Li is the free load and 
conservation gives us the bound load as Lit— Lj. 
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These four equations are presented in de-dimensionahzed form, 
with u,v,£i,£2 representing the dimensionless concentrations of 
Repressor 1, Repressor 2, Loadl and Load2 respectively and x the 
de-dimensionalized time. The basal parameter values that we use 
are as follows: a; = otj = 0.2; fS/ = fij' = 4; n = 3; k^ni' = k„n2' = 0.5; 
koin' =k„fi2' = 0.5;ki=k2=l; [Lit] and [Ljt] are variable. Note 
that Equations (1) and (2) without the last two terms incorporating 
the load are the standard equations for analyzing the toggle switch 
that have been widely used in both empirical and theoretical work 
[5,36]. These equations are discussed in more detail in Supple- 
mentary Text S 1 Section 1.1. The derivation of this model follows 
that of Kobayashi et al [33]. AH parameters excluding load binding 
rates were sourced from Kobayashi et al [36] ; extensive parameter 
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Figure 1 . Schematic diagram of the circuits studied in this paper. (A). The basic toggle switch is the network shown without the dotted line. 
Repressor 1 represses the production of Repressor 2 and vice versa. The dotted line denotes a positive feedback motif found in some natural circuits. 
(B). A cartoon of part of the IVIAPK activation pathway in T lymphocytes, adapted from [29], showing the role of Ras activation. Signals from peptide- 
IVIHC complexes are received at the TCR and lead to phosphorylation of the cytoplasmic chains of the TCR by the Src kinase, Lck. This recruits the 
kinase ZAP70 which trans-autoactivates and phosphorylates a scaffold called LAT, which recruits Grb2 and SOS to the plasma membrane. SOS 
activates Ras as as shown. (C) A simplified model of the Ras switch. RasGDP transforms into RasGTP via the enzyme SOS. However the catalytic rate of 
SOS increases when bound to RasGTP. This sets up an autocatalytic positive feedback. RasGTP is deactivated by enzymes called RasGAP's (among 
others). 

doi:1 0.1 371/journal.pcbi.1 003533.g001 



sensitivity of the load binding rates was performed and are discussed 
in the Supplementary Text S 1 section 1 .4 and Figs. S 1 , S2, Table S 1 
and Figs. S15 and SI 6. The effect of a load arises from tlie binding 
competition between the promoter wliere the repressor binds and 
the load. This competition is not directly incorporated into tlie HUl 
function, since the binding step witli the promoter is not explicitly 
modeled and is treated in an effective way. In reality however the 
concentration of the promoter is so small compared to tliat of the 
load, that the use of HiU functions is justifiable [37]. There are 
possibly exceptional cases such as a high copy number of plasmids 
compared to load concentrations where this assumption does not 
apply. Note that the HiU function is an effective phenomenological 
equation describing gene transcription and protein production, and 
standard Law of Mass Action (LMA) methods to derive the HiU 
functional form may not apply for many transcription factors that 
nevertheless show HiU kinetics [38] . Thus it is preferable to use HiU 
function forms for this analysis. 

To calculate transition times, we first start the system in one 
state, say high Repressor 1. After the system has reached steady- 
state, we add a constant concentration of the inducer and measure 
the time taken for Repressor 2 to go from 10% of its maximum 
value to 90% of its maximum value. This is the "rise time". 
Similarly the "decay time" is the time taken for Repressor 1 to go 
from 90% of its maximum value to 10% of its maximum value. 
The level of the inducer remains fuced. 

In practice the inducer may decay and the transition would 
depend upon there being inducer present for a sufficiently long 



time to induce transition. In such cases the amount of inducer 
required may be of interest. When the inducer is apphed as a bolus 
with a first order decay rate, it appears as an exponentially 
decaying pulse. We thus included a fifth differential equation 
governing the amount of Inducer. 



Here dj is the ratio of the inducer degradation constant to the 
repressor degradation constant. We used Eq. 5 only when 
estimating the amount of inducer required to switch states for 
different loads and different decay rates of the load (Supplemen- 
tary Text SI section 1.4 and Supplementary Tables S3, S4). 

A genetic toggle switch can be induced to change states by the 
alternative method of repressing the highly expressed repressor, 
and in fact the original toggle switch used this form of induction 
[5,33]. We repeated our calculations for the basic model for the 
case of alternative induction, but found no qualitative differences. 
The alternative induction model along with the equations is 
detailed in the Supplementary Text SI section 1.4. 

Equations 1-4 assume that the load itself stays in steady state 
during the switching of the toggle between one state and another. 
However in reality if the load is another protein, it is also 
synthetized and degraded by the cell, and therefore its level could 
be dynamic. We also simulated this situation by incorporating a 
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synthesis and a degradation rate for each load. This resulted in 
Equations 3 and 4 being replaced by: 



-r- = -koni ki —u£i+k„ffi ki- — (ci)+ — —£i (6) 

ai kd2 kd\ o 0 



Ml , kbi I, , , kb2 1 s kd2 kd2 1, 

= -ko„2 K2 7 Vil+k„ff2 k2-, — (C2) + ^ ^ ll (7) 



dx 



Here c\ is the load-repressor complex and k^i and kd\ are the 
synthesis and degradation rates respectively for Repressor 1, and 
correspondingly for Repressor 2. The parameters are defined in 
the Supplementary Text SI, section 1.5. Since the total load is no 
longer conserved, we need to include additional equations for the 
load repressor complex. 



dc\ 
dx 

dC2 

dx 



= ko„i kiuii —koffi k\ci 



= Knl k2Ul2 - koff2 k2C2 



(8) 



(9) 



Our model assumes that when the repressor protein is bound to 
the load, it is protected from degradation. However it is also 
possible that even when the protein is bound to the load, it can still 
degrade. To check the impact of removing the protection 

assumption, we also consider an additional model where the 
repressor can still degrade with the same rate constant when 
bound to the load. The equations for that model are shghdy 
modified versions of the equation above, and are presented in 
detail in the Supplementary Text SI, section 3.2. 

We conducted parameter sensitivit)' analysis on models utilizing 
both forms of induction; these did not show any qualitative change 
on wide variation of key parameters (Tables SI, S2, S3, S4 and 
Supplementary Text SI). 

Toggle with positive feedback 

A positive feedback was added to the Rl side of the toggle 
switch as an inducible promoter with a Hill coelficient of 1. We 
assumed that the positive feedback acted on the same promoter as 
the repression, resulting in a composite term for production of Rl 
from promoter 1 where p is the strength of positive feedback. 



dRi 
dt 



\+Ri/k5 + FC^/k"i 



-diRi 



(10) 



The derivation of this equation can be found in the 
Supplementary Text SI, section 1,6.1. As before, ai is the leaky 
production of Ri while oti+Pi represents the activity of the 
promoter in the absence of repression or positive feedback. We 
chose k2 and k-, = 1 , di = 0.2, and for the figures in the main paper 
we chose p = 3.5. We address other values of the positive feedback 
in Fig. S6 and the Supplementary Text SI, section 1.6.2. 

Stochastic simulations 

We perform stochastic simulations and histogram the concen- 
trations of the repressor proteins to construct their probability 
distribution. The quasi-potential of the toggle is given by the 



negative logarithm of this probability distribution [39]. In order to 
construct the probability distribution we make use of the 
phenomenon of noise-induced switching. Recent theoretical work 

has shown that multiplicative noises due to stochastic fluctuations 
can induce switching [40-42]. Experimental results demonstrate 
bimodal populations that correspond with theoretic predictions 
arising from noise-induced switching [41]. 

Stochastic simulations were carried out using a modified Gillespie 
algorithm using the standard rate expressions for every reaction 
(Table S5). We chose a reaction volume that would correspond to a 
small number of molecules in the system. Stochastic fluctuations 
then drive the system to transition between states rapidly, allowing 
us to collect sufiicient data points. In order to make sure that the 
system was not being biased by the small volume, we also repeated 
the calculations for a five times larger volume (and hence molecule 
number) and found quahtatively similar results (Fig. S4). 

For the positive feedback toggle switch the same equations were 
used except for the repressible production of Repressor 1, where 
we used instead the rate expression gi\ en by the right hand side of 
Eq. 10 in the Monte Carlo simulations. 

Ras-kinase system 

For our study we adapted the minimal model of the Ras switch 
proposed by Das et. al. [30] with the addition of a reversibly 
binding load in the form of the Raf protein (Fig. IC). The model 
contains three proteins, Ras, which exists as RasGDP or RasGTP, 
SOS, the guanine exchange factor (GEF) that catalyzes the 
transformation from RasGDP to RasGTP and a GTPase, 
RasGAP. SOS on its own has very low GEF activities. However, 
the activity of the GEF pocket is strongly influenced by the binding 
state of an allosteric pocket in Cdc25 domain [29,30]. When the 
allosteric pocket is bound by RasGDP, the GEF activity is 
increased by 5 times. If the allosteric pocket is bound by RasGTP, 
its GEF activity is increased by 75 times. In this way, RasGTP can 
upregulate its own production rate by binding to SOS, thus 
constituting a positive feedback loop. RasGTP is deactivated by 
GTPase's such as RasGAPs that are constitutively present. 

After Raf binds RasGTP, the complex catalyzes the phosphor- 
ylation of Raf leading to a phosphorylation cascade. For this study 
we ignore Raf activation and only consider the effects of Raf as a 
binding partner for RasGTP. The Das paper [30] also models the 
systems using Michaelis-Menten (MM) forms for the actions of the 
enzymes which is quite standard for modeling systems of enzymatic 
reactions. However since in this model the load competes not with a 
promoter, as in the toggle switch, but with another protein, it is 
possible that the quasi-steady state assumption of the MM form 
could be introducing some inaccuracies in the results. To account 
for this possibility we wrote the entire model using the Law of Mass 
Action. We separately simulated the model using the MM func- 
tional forms (Supplementary Text S 1 section 2 and Figs. S7 and S9). 
The equations for the MM forms are listed and discussed in detail in 
the Supplementary Text SI. The reactions and rate constants for 
this model are listed in Table S6 and Table S7. 

We use the following notations for the species involved in the 
system: 

xi = [SOScat] ; X2 = [RasGDP] ; Xi = [RasGTP] ; 

X4 = [SOScat{RasGDP)] ; jcs = [SOScat{RasGTP)] ; 

xe = [SOScat{RasGDP) : RasGDP]; 

X-, = [SOScat{RasGTP) : RasGDP] ; = [RasGAP] ; 

X9 = [RasGAP : RasGTP] ; xio = [Raf] ; xn = [RasGTP : Raf] 
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dt 
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Moreover, four of the basic protein species along with the 
complexes they participate in have associated conservation laws. 
These are as follows: 



SOST = x\+X4+Xi+Xd+x^ 



(22) 



RasT = X2+x-i+X4+xs + 2x6 + 2x7 + X9 + xi 1 (23) 



GAPT = Xi+X9 
J?a/r=xio + xii 



(24) 
(25) 



In the Ras model too we implicitly assume that when RasGTP 
is bound to Raf, it is protected from de-activation by a RasGAP. 



We also study the effects of relaxing this assumption on both the 
LMA and the PSSA models. The modifications to the original 
model are detailed in the Supplementary Text SI section 3.3.1. 

We used XPPaut to perform a bifurcation analysis of the Ras 
switch with changing levels of SOS, with and without a load. The 
quasi-potential landscape does not provide useful insights into load 
induced modulation of the Ras switch and hence the probability 
distributions are not reported. 

Results 

The bistability properties of the toggle switch do not 
change unless the repressor can degrade when bound to 
the load 

The presence of a binding partner for either Repressor 1 or 
Repressor 2 (which we refer to thereafter as the load) introduces 
new terms in the differential equations describing the toggle 
switch, i.e. the last two terms in Eq. 1 and in Eq. 2, as well as two 
new equations, Eq. 3 and 4, in the dynamical system. However it 
can be easily seen that in steady state Eq. 3 and 4 are also 
independently set to zero, and therefore do not affect the 
bifurcation properties of the switch. Even in the case of a dynamic 
load, since Eq. S13 and S14 are set to zero to ensure the load- 
repressor complex is in steady state, the additional terms in Eq. S9 
and SIO are also zero. Thus the load makes no difference to either 
the bistability of the switch or to the parameter values where the 
bistability is seen. 

The exception is when the repressor molecule can degrade even 
when bound to the load, which may be relevant in some 
experimental situations. As Fig. 2A shows, when a load is added 
symmetrically to both sides of the toggle switch, the two stable 
states approach each other and eventually annihilate, leaving a 
monostable system. Fig. 2B shows that when a load is added only 
to one side, the system again goes from bistable to monostable at 
some critical value of the load. In effect, the upper stable point 
vanishes and is no longer accessible due to leakage of the repressor 
affected by the load. 

The reason for the change in steady state behavior is made clear 
on examining the equations of the system. Here we need to 
incorporate additional reactions that represent the decay of the 
repressor-load complex into the load alone. This k'ads to an 
additional term in the equation for the load and the repressor-load 
complex (Eq. 844 and S45). However this term does not appear in 
the equation for the repressors, which continue to be governed by 
Eq. 1 and Eq. 2. As a consequence in the steady state, the 
additional terms in Ecj. 1 and 2 no longer equal zero and the 
steady state properties of the switch are influenced by the presence 
of the load. 

As can be seen from an examination of the chemical reaction 
system, this mechanism of abrogation of bistability arises whenever 
the load-repressor complex participates in a non-reversible (from 
the repressor's point of view) chemical process that leads to an 
unbalanced leakage of the repressor from its function as a 
repressor by the presence of the load. A more interesting example 
of such a process could be provided by a chemical reaction system 
where the load is an enzyme for one of the repressor molecules, 
which is transformed by the enzymatic action into a protein no 
longer capable of repression. The mathematical analysis of this 
case is exactly the same as the model we are currendy discussing 
hence we do not consider it separately here. 

However a load can significantly change the dynamic response 
of the basic genetic toggle switch as we shall see below. We exam- 
ined two different measures of dynamic response, response time 
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Figure 2. Bifurcation diagram of the genetic toggle switcii wKien the repressor can decay from the ioad-repressor complex. The thick 
lines are stable steady states, the dashed lines are unstable steady states. (A). A load is added symmetrically to both sides of the toggle. The stable 
states of only one Repressor molecule with respect to the load are shown. With zero load the toggle switch is bistable with well separated steady 
states. As the load increases, the two stable states approach each other and the unstable state, and eventually merge in a bifurcation at a critical 
value of the load. The system is monostable beyond this critical value. (B) A load is added only to Repressor 1. The high state of Repressor 1 
approaches the unstable steady state as the load increases and merges with it at a critical value of the load, leaving only the lower state accessible to 
the system. 

doi:10.1371/journal.pcbi.1003533.g002 



for state switching and the amount of inducer required for .state 
switching. 

The response time for state switching of the toggle 
switch increases 

We measured two response times, the rise time which quantifies 
the time taken for the concentration of Repressor 2 to increase 
from its low or zero level in state 1 to its high level in state 2, and 
the decay time which measures the time taken for Repressor 1 to 
decay from its high level in state 1 to its low level in state 2, in both 
cases in response to a constant inducer. Specifically the rise time 
measured the time to go from 10% to 90% of the steady state 
maximum, while the decay time measured the time to go from 



90% to 10% of the steady state maximum. These measurements 
were made using the deterministic model in the cases when the 
load was applied only to one side and to both sides of the switch. 

We found that both the rise time and the decay time increase 
with increasing load concentration. Interestingly, this relationship 
was approximately linear in all cases (Fig. 3 A & B). The slope of 
the linear relationship represents the increase in response time due 
to unit increase in load. We found that the slope of the line was 
larger when the load was applied to the opposite side of the system 
before the switching rather than the same side (Fig. 3A), indicating 
that it is harder to switch out of a state without a load to a state 
with a load than the reverse. However when a load was applied to 
both sides, the slope of the linear fit was higher than when the load 




Figure 3. Effects of a load on transition times of the basic toggle switch. (A). The time taken to reach 90% of maximum value for the protein 
undergoing a low-to-high transition as a function of the Load, normalized by the steady-state amount of Repressor 1. Normalized time is a unit-less 
number defined by the transition time (rise or decay) of the system at a given loading condition divided by the transition time (rise or decay) of an 
unloaded system. (B). The time taken for the concentration of the protein undergoing a high-to-low transition to reach 10% of its maximum value. 
The X- and y- axes are the same as for the previous panel. 
doi:1 0.1 371/journal.pcbi.1 003533.g003 
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was only on the opposite side, suggesting that both the "opposite 
side" and the "same side" delays are operating. 

While we also found an approximately linear relationship 
between the decay time and the concentration of the load, there 
was little difference between the decay times for the state with the 
load ("same side load") and the state without a load ("opposite side 
load") at our base parameter values. Thus the load affects rise time 
and decay time differently. When a load was applied to both sides 
of the switch, the slope of the decay time linear fit was larger, again 
indicating the operation of both delays. 

We tested these results by changing parameter values for the 
binding of the load (Table SI) and found that in all cases we obtain 
a good linear fit for the response time. For the rise time, the slope 
was uniformly larger when the load was applied to the opposite 
side as compared to the same side, and it was the largest when 
loads were applied on both side.s. For the decay time, the slope 
could be larger or smaller when the load was applied to the 
opposite side of the decaying state compared with the same side, 
but it was always larger than both when a load was applied on 
both sides. The slope depended non-monotonicaUy upon the 
dissociation constant (Kd) of the binding between the repressor 
protein and the load, with both low Kd and high Kd having a 
smaller effect that those in between (Fig. SI). This was because 
when the Kd was low, i.e. strong binding, the concentration of the 
load-repressor complex was unaffected by the state of the switch. 
However when the Kd was high, the maximum concentration of 
the load-repressor complex was smaller, thereby having a lesser 
effect on the system (Fig. S2). Thus response times are maximized 
when the load acts as a dynamic sink, i.e. it takes up newly 
synthesized repressor when the state changes from the unloaded to 
the loaded side, and releases the bound repressor when switching 
from the loaded side to the unloaded side. 

Previous studies of response times of biochemical networks with 
and without a load have also seen monotonic increases in the 
response time of simple transcriptional circuits [37]. However the 
extremely consistent approximately linear response we see under 
wide variation in parameter values is extremely intriguing. 

An increase in response time should also imply that the 
concentration of inducer required to shift states should also be 
affected, especially when it can decay. In atxxjrdance with this 
expectation we also found that the concentration of inducer 
required to switch states increased exponentially with increasing 
load, as seen in Table S2. The parameter of the exponential fit was 
dependent on the inducer decay rate, indicating that the amount 
of time the inducer remains abcn e a threshold is the key factor 
governing the switching. We find that this response to a load is 
unaffected by the mode of switching the toggle, and induction by 
repression of the current state yields the same qualitative results 
(Table S2 & S3). 

In our analysis so far we have assumed that the total 
concentration of the load is fixed. We now analyze the case when 
the load is generated by a constitutively active promoter and can 
decay at a first order rate. We find that in this case too the 
qualitative features of the transition time remain the same as the 
toggle switch with a fixed load, i.e. it was approximately linear in 
all parameter regimes tested (Supplementary Text SI section 1.5, 
Fig. S3 and Table S4). The reason why we do not see a difference 
from the basic toggle switch is that the transition times ultimately 
measures time between steady states, and we wait for the system to 
come quite close to the steady state value (90%). Thus the 
concentration of the load has also reached a steady state value and 
the system behaves as it would with a fixed load. 

We also tested the response times when the repressor can leak 
away from the systems after binding with the load. Here we find 



that (Fig. 4) when a load is apphed to the same side, the rise time 
continues to increase monotonicaUy linearly with the load but the 
decay times decreases monotonicaUy with the load. However when 
a load is applied to both sides, we find a negative linear relation 
between the transition times for both rise and decay and the load. 

The reasons for the change in behavior is because as we saw 
previously, when the repressor can leak away from the repressor- 
load complex, a load has a dramatic effect on the bistabiHty 
properties of the switch, abrogating bistabUity very quickly (Fig. 1). 
When only one repressor has a load, the high state of that 
repressor approaches the unstable state, indicating a decrease in 
the domain of attraction. Shifting out of that state thus becomes 
easier with increasing load. When both sides have loads, both 
stable states approach the unstable state, therefore shifting out of 
either state becomes easier, and both transition times decrease. 

Dramatic changes in the potential energy landscape and 
probability distributions of the toggle switch 

The modulation in the dynamic properties of the basic genetic 
toggle switch discuss(;d ab()\'e suggests that the load has altered the 
potential energy landscape of the toggle switch, making it harder 
to switch. For two-dimensional and higher systems, such as the 
toggle switch, analytical methods to construct the potential 
landscape are not available, but a quasi-potential can be 
constructed from the probability distribution function of the 
concentrations of the repressor molecules, where the quasi- 
potential is given by the negative of the natural logarithm of the 
probability distribution [43,44]. To calculate this we performed 
Monte Carlo simulations of the toggle switch using a Gillespie type 
algorithm elaborated in the Methods section. When the toggle 
switch is symmetrically balanced, both the probability distribution 
function and the potential energy landscape are completely 
symmetric. If the system is started in State 1, random fluctuations 
can drive it into State 2 and vice versa. The probability 
distribution can then be constructed by counting the frequencies 
of these random fluctuations. However since the genetic toggle 
switch can be very stable, a numerical computation of the potential 
energy landscape requires impractically long simulation times (as 
we show below). While computational methods to sample rare 
trajectories in such cases exist, they are very sensitive to choices of 
parameters [42,45]. We developed a computational protocol in 
order to numerically obtain the probability distribution function of 
both protein concentrations and the transition times. We chose an 
appropriate volume for the genetic toggle switch such that exactiy 
the same parameters as in the deterministic simulations led to the 
operation of the toggle switch with only a small number of 
proteins. The toggle remains bistable in this regime but the small 
protein numbers vastiy increases spontaneous stochastic fluctua- 
tions arising out of multiplicative noise in the system and allows the 
simulation to explore parameter space and collect enough data. 

Our simulations showed that the switch switched states a large 
number of times. In order to account for diflerences in the time 
step in diflerent states, the probability density function of the 
concentrations was constructed using a time trace collected after 
approximately 1 second intervals. As Fig. 5 shows, for a symmetric 
switch we obtain a symmetric bi-modal probability distribution 
that corresponds to a double-well potential. 

When we add a load to the system asymmetrically, in the form 
of a binding partner for the Repressor 1, we find that the 
probability distribution becomes extremely skewed, and the total 
weight of the probability distribution corresponding to the other 
side, i.e. Repressor 2, dramatically increases (Fig. 6A). This 
indicates that the underlying double well potential has become 
skewed and the state 2, corresponding to high Repressor 2, has 
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Figure 4. Effects of a load on transition times of a toggle switch without the protection assumption. (A). The time taken to reach 90% of 
maximum value for the protein undergoing a low-to-high transition as a function of the load. The system is de-dimensionalized as described in 
Supplementary Text SI section 1.1 and 3.2.1. (B). The time taken for the concentration of the protein undergoing a high-to-low transition to reach 
1 0% of its maximum value. Note that the linear relationship for both-sided load transition times, and same-sided decay time, and opposite-sided rise 
time has a negative slope. The relationship for same-sided rise time and opposite sided decay time has a very small, but positive slope. 
doi:10.1371/journal.pcbi.1003533.g004 



increased its stability at the cost of State 1 (Fig. 6C). When a load is 
apphed to both sides symmetrically, the concentration probability 
distribution reverts to a symmetric bimodal distribution corre- 
sponding to a symmetric double-well potential (Fig. 6B & D). 

In order to test this directly we calculated the distribution of 
lifetimes in state 1 and the lifetimes in state 2. As shown in Fig. 7, 
when the switch is symmetric with no load, the lifetime distribution 
is exponential, as should be expected for a simple two-state system. 
However when the load is applied to Repressor 1 , the probability 
distribution of the lifetime in state 2 increases dramatically. The 
average lifetime of state 1 also increases but only by a very small 
amount. The time spent in state 2 does not appear to saturate, and 
continues to increase with increasing load. When loads are applied 



symmetrically to both sides, the lifetime histogram in Fig. 7 indi- 
cates that both sides have been stabilized since the system spends 
significandy longer time in each state. Note that in an equilibrium 
system this would have been indicated by the deepening of the 
potential well. However in non-equilibrium systems the potential 
well picture does not completely capture the dynamics and there is 
an additional contribution from a "curl flux" [43,46] that needs to 
be taken into account. For our purposes calculating both the 
distribution of concentrations and the distributions of lifetimes 
captures the dynamics of the toggle switch. 

To test whether our results change for higher protein con- 
centrations, we increased protein concentrations about fivefold 
and recalculated the probability distribution function. We find that 




Figure 5. The probability distribution function and the quasi-potential of the genetic toggle switch without a load. (A). The 
probability distribution function of a toggle switch without a load. The x- and y- axes here represent the number of molecules of Repressor 1 and 
Repressor 2 respectively, while the z-axis is the frequency of its occurrence. Note that the distribution is symmetric as expected. (B). The quasi- 
potential of the symmetric toggle switch, showing the symmetric double-well potential constructed by taking the negative logarithm of the 
probabilities shown in (A). A small offset of 0.001 was added to the probabilities to prevent taking the logarithm of zero. This does not change the 
shape of the well. 

doi:1 0.1 371/journal.pcbi.1 003533.g005 
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Figure 6. The probability distribution function and quasi-potential of a toggle switch with a load. The 3-dimensional plot is viewed with 
the xy-plane horizontal for better contrast. The x- and y-axes are numbers of molecules of R1 and R2 while the z-axis is either probabilities or the 
quasi-potential. (A). The probability distribution function (pdf) of the toggle switch of Fig. 5 but now with a load of 20 molecules on Repressor 1 (R1 ). 
(B). The pdf of the toggle switch with a load of 20 molecules on R1 and 20 molecules on R2. (C). The quasi-potential of the toggle with a load of 20 
molecules on R1, i.e. corresponding to panel A. (D). The quasi-potential of the toggle with equal loads of 20 molecules on each repressor, i.e. 
corresponding to panel B. 
doi:10.1371/journal.pcbi.1003533.g006 



our qualitative results remain robust despite the increase in protein 
concentrations (Supplementary Text SI section 1.3 and Fig. S4). 
Switching between states is rare at these protein numbers, with a 
mean residence time in state Rl for the unloaded switch being 
approximately 6x lO' min against about 700 min for the basal 
case considered, a diflFerence of almost three orders of magnitude. 
However as for the basal case, the quasi-potential landscape skews 
significantiy with the addition of a load on the switch. 

"Opposite Side effect" dominates the load effect in the 
basic toggle switch 

These results allow us to interpret the dynamic results that we 
obtained earlier. If the system is in state 2 and there is a load on 
state 1, a transition requires an increase in Repressor 1 concen- 
tration in order to suppress the production of Repressor 2. A load 
on Repressor 1 however competes with the promoter of Repressor 
2 for binding with Repressor 1 , and thereby reduces the effective 
concentration of Repressor 1. This effectively stabilizes state 2. 
The dynamic analysis shows that state 1 not only remains an 
attractor state but in fact it takes a longer time, and more inducer, 
to shift out of state 1 as compared with the no-load situation. This 
is because the load also acts as a reservoir for Repressor 1 , and in 
fact increases its total concentration. This slows down the tran- 
sition to state 2. Interestingly this "same side effect" is generally 
weaker than the "opposite side effect" above. In agreement with 
this picture, the stochastic simulations show that the distributions 
of lifetimes in state 1 broaden slighdy on addition of a load. 

If the load is present symmetrically on both sides, the concen- 
tration histograms in Fig. 6 and the time histograms in Fig. 7 indi- 
cate that both states have been stabilized, due to a combination of 
the 'same side' and the 'opposite side' effect now acting together to 
stabilize each state of the switch. In the dynamical simulations this 



is seen by the increased slope of the response time line for the case 
of a load on both sides. Results for additional parameter values are 
shown in Fig S15 and Fig SI 6. 

Positive feedback moiety makes toggle switch tunable 

When a positive feedback moiety is introduced in the toggle 
switch, we again see a linear relationship between the rise time and 
the decay time of the two states of the switch and the load (Fig. S5). 
Therefore here too the load appears to be skewing the underlying 
potential landscape of the switch. Using stochastic simulations we 
constructed the probability distribution function of this toggle 
switch as described above. We found that even in the absence of a 
load, when a positive feedback moiety is introduced on one side of 
a toggle switch, the probability distribution for the toggle switch, 
and hence the quasi-potential landscape, becomes extremely 
skewed in favor of the state with positive feedback as shown in 
Fig. 8A. Even with no load on the system, the switch is biased to 
State 1 and the lifetime spent in State 1 is much longer than in 
State 2. If a load is added to R2, the opposite side effect 
additionally favors State 1. If a load is added to Rl however, the 
opposite side effect favors State 2 (Fig. 8B). It is possible to balance 
these effects resulting in a more even distribution by adjusting the 
load on Rl and the strength of positive feedback. As the load on 
Rl is increased beyond this balance point, the opposite side effect 
dominates and the probability distribution becomes skewed 
toward State 2 (Fig. 8C). As the opposite side effect increases 
with increasing load, the lifetime in State 2 also increases in 
agreement with the findings for the regular toggle switch (Fig. 8D). 
The lifetime in State 1 also increases by a smaller amount, as for 
the regular toggle switch (Fig. 8E). 

For the toggle switch with the positive feedback moiety, we can 
also check the consequences of allowing repressor leakage through 
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Figure 7. Distribution of the lifetimes of the toggle switch with and without loads. The time the system spent in either state R1 or state R2 
was calculated from the time trace of the stochastic simulations and a histogram made of the results. The histogram is shown on a semi-log plot to 
accommodate the data on a single chart. (A). Lifetimes in State R1. The unloaded state is the solid curve that is to the extreme left of the others, 
showing that the lifetimes in state R1 increase slightly on addition of load on R1 alone due to the "same side effect". (B) Lifetimes in State R2 when 
load is on R1. The solid curve on the extreme left is the unloaded state. There is a significant increase in lifetimes due to the "opposite-side effect" of 
the load on R1. (C). Lifetimes with a balanced load, showing that both the states R1 and R2 get stabilized with a significant increase in lifetimes on 
addition of a small load on both sides. Note that the distributions for Rl and R2 for equivalent cases coincide as should be expected. 
doi:1 0.1 371 /journal.pcbi.1 003533.g007 



the repressor-load complex. As shown in Fig. SI 3, this addition to 
the system affects the steady state properties of the switch and 
bistability is abrogated after the load increases beyond a critical 
value, when load is present for both sides or only one side. 

Loads fundamentally transform positive feedback based 
switches in signal transduction 

The RasGTP system shows a bistable transition from a low 
RasGTP state to a high RasGTP state as the activating signal, in 
our case the number of SOS molecules, are varied. As Fig. 9 
shows, a system with no Raf shows a classic Z-shaped bifurcation 
diagram with two bifurcations as SOS is varied. The first bifur- 
cation marks the transition from a monostable low-RasGTP state 
to a bistable system with a "high" RasGTP state (and an unstable 
intermediate state). The second bifurcation marks the transition 
from the bistable state to another monostable state with a high 
concentration of RasGTP. 

When Raf is added to the system, the bifurcation diagram 
changes and the two bifurcations start approaching each other. 
This is because the effect of adding Raf is equivalent to seques- 
tering away some of the activated RasGTP in an "inactive" com- 
plex. When Raf concentration crosses a threshold, the bifurcations 
annihilate each other and disappear. This system is now char- 
acterized by a single stable point for all concentrations of SOS, 
and the disappearance of the threshold for Ras activation. While 



there appears very litde free Ras, in reality, even for low SOS 
concentrations there is a large concentration of the activated 
RasGTP-Raf complex (since RasGTP in these complexes is also 
protected from the action of the Ras GTPases). 

This can be seen in another way in Fig. S8 where the stable state 
of RasGTP is plotted against the level of total Raf in the system, 
keeping the level of SOS constant. Again we see that a bistable 
system is transformed into a monostable system when Raf increases 
beyond a threshold. These results are exactly the same for the model 
which assumes Michaelis-Menten kinetics, except for small changes 
in molecule numbers, as can be seen in Fig. S7 and 89. Results do 
not change on changing load-binding parameters (Fig. SIO, SI 1) 

Thus the addition of the Raf scaffold, which is an integral part of 
the MAPK cascade, fundamentally changes the qualitative behavior 
of the positive feedback switch. The main reason why the steady 
state bifurcation properties are affected here in contrast to the basic 
genetic toggle switch is that for this signaling circuit, as seen in Eq. 
22-25, total Raf and Ras are conserved, as is typical for a short 
timescale signal transduction system. These conservation laws 
couple Raf concentration to RasGTP concentration even at steady 
state. Therefore adding Raf to the system effectively reduces total 
Ras concentration since Raf sequesters away Ras from the switch. 

To see this more generally, consider for example a chemical 
reaction system comprising of n-species Yi,...Y„. Let us assume 
without loss of generahty that the species ¥„ is coupled to a 
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Figure 8. The genetic toggle switch with a positive feedbacl< motif on Repressor 1 (R1 ). (A). The probability distribution function (pdf) with 
no load. The positive feedback on Repressor 1 leads to a pdf skewed in favor of R1. (B). The pdf with a load of 20 molecules on R1 showing the 
increase in the weight of R2 due to the "opposite side effect". (C). The pdf with a load of 40 molecules on R1. This load is more than enough to skew 
the pdf in favor of state R2. (D). Histogram of lifetimes in R1 with varying levels of load on R1 . Comparison with panel A shows that the unloaded state 
has been stabilized by the positive feedback. Note that the lifetimes increase very slightly due to the "same side effect". (E). Histogram of lifetimes in 
R2 with varying levels of load on R1 . The unloaded case is the curve on the extreme left. Note the initial asymmetry in the lifetime distribution due to 
the positive feedback, as well as the large increase in lifetimes with the inclusion of a load. 
doi:10.1371/journal.pcbi.1003533.g008 



downstream circuit through a binding reaction with a load, L. The 
(n+2) differential equations describing this system are: 



[YnL] 



(30) 



dt 



--fi(Y,,...Y„) 



(26) 



dY„ 



--f„(Yi ,...,Y„) + k„„[Y„L]- Kff [ r„] [L] (27) 
d[L\ 



dt 



--koff[Y„L]-k„„[Y„][L] 



d[Y„L] 
dt 



--koff[Y„L]-ko„[Y„][L] 



(28) 



(29) 



Note that for simplicity of notation we have not indicated the 
dependence of the dynamical system on its own parameter values. 
Now in the steady state, if the set of equations is complete, the left 
side uniformly goes to zero and we recover the result that the 
steady state remains exactly the same with or without a load, as for 
the genetic toggle switch. However let us now assume that we have 
an additional conservation law, say. 



This conservation law implies that one equation in our 
dynamical system is redundant, and we need to drop one equation 
to make the system linearly independent. We can decide to drop 
Eq. 19, and substitute Y„=Y°-[Y„L] in Eq. 20 and Eq. 21 and 
solve the resulting (n+1) equations for the (n+1) unknowns, 
Y],...,Y„_l,Y„L,L, obtaining Y„ as a residual from Eq. 22. Thus 
the steady state solutions of the F/s now involve the amount of the 
load. Clearly, the existence of the conservation law has led to a 
change in the steady state properties of the dynamical system. 
Note that Y„ itself would usually enter (by itself or in the form of 
other complexes, which then would also need to be accounted for in 
the conservation law Eq. 22) into one or more of the equations for 
the remaining species, Fi,...F„_i. This would result in the 
equations for those other species explicitly involving, and thus 
depending upon the level of the load. For the Ras system above, Eq. 
16 couples the load, Raf, to the concentration of Ras. However Ras 
concentration and SOS concentration are also coupled. Thus the 
load explicitly affects the steady state values of all species 
concentrations in this system. This leads to a fundamental qua- 
litative change in the bifurcation properties of the system. 

Discussion 

It has been pointed out previously that significant sequestration 
effects can abrogate zero order ultrasensitivity [26,47,48], can 
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Figure 9. Bifurcation diagram of the Ras switch with different 
levels of Raf (load) on the system. The total number of SOS in the 
simulation box is used as the parameter being tuned, which varies from 
0 to 1 000. For Raf = 0, Raf = 1 0 and Raf = 30, there are two bifurcations as 
SOS is increased. In the first bifurcation a new high valued stable steady 
state appears along with the low valued stable steady state. In the 
second bifurcation, the low valued stable state disappears leaving 
behind only the high valued state. The dotted line marks the unstable 
steady state that also comes into existence in the bistable region. As 
total Raf increases, the two bifurcations approach each other. When 
Raf =50, the system has lost both of its bifurcations and is characterized 
by a single stable steady state at all values of Raf. 
doi:1 0.1 371/journal.pcbi.1 003533.g009 

change the dynamics of simple phosphorylation circuits [23,24] 
and change oscillatory behavior in some circuits [27]. We add to 
this body of work by demonstrating that the addition of a simple 
binding partner to the output protein of a genetic or signaling 
switch can have dramatic effects on its properties, and can 
fundamentally change the operation of the switch. 

For a genetic toggle switch with two mutually repressing 
proteins such as the classic switch built by Gardner et al. [5] we 
showed that even though the presence of the binding partner does 
not alter steady state properties of the switch, it can drastically 
change the dynamic properties. Using a novel potential landscape 
analysis, we showed that this is because the addition of the binding 
partner skews the underlying quasi-potential, making one state 
significantly more stable than the other. In practice therefore, a 
genetic toggle switch that is significantly skewed towards one side 
may never properly function as a switch. Thus the downstream 
consequences of such loads need to be taken into account when 
designing larger synthetic circuits with the toggle switch as one of 
the elements. 

On the other hand this phenomenon actually provides a way of 
making artificial switches tunable. It is possible to engineer a 
biased switch merely by adding a load on the opposite side of the 
toggle, which is a useful device when engineering a switch that is 
designed to be switched on only in special circumstances. A load 
on both repressor proteins similarly stabilizes both sides of the 
toggle switch. This could be useful when working with synthetic 
components with low concentrations in cells, especially those that 
display stochastic switching. A load on both repressor proteins can 
significantly increase the stability of such a toggle. 



In natural systems, mutually repressing toggle switches are often 
found with other complexities, such as a positive feedback motif on 
one side. The positive feedback motif by itself biases the toggle 
switch by stabilizing the side it is on at the expense of the other 
side. A load on the same side then stabihzes the opposite side, and 
can re-establish balance between the two quasi-potential weUs. For 
engineering circuits in multi-cellular organisms, it is worth noting 
that that feedback between the load on a toggle switch and the 
strength of the positive feedback may ensure that the switch 
operates efficiently even in the presence of cell to cell variability in 
the load. How loads vary between cells and in multi-cellular 
organisms is an interesting question to explore in future work. The 
presence of the positive feedback provides a potential target for 
evolutionary fine-tuning of the switch. 

In the above analyses we use novel potential landscape methods 
that have proved useful and insightful in fields such as protein 
folding to discuss the fundamental properties of a dynamical 
system that shows not apparent changes in its stability properties. 
We demonstrate that these methods, though still relatively 
underdeveloped for use with non-equilibrium chemical reaction 
systems, hold promise for understanding the dynamics of such 
systems beyond what linear stability analysis can provide. How- 
ever there are certain conditions when addition of a load changes 
the stability properties of the genetic toggle switch. One class of 
such effects happen when the repressor can leak away from the 
repressor-load complex, as can happen either when the repressor 
can decay or degrade when bound to the load, or when the load 
can modify the repressor and make it unable to repress. We show, 
employing standard bifurcation analysis, that additional loads in 
this system can abrogate the switch-like properties of the toggle 
switch entirely. 

In switches based on autocatalysis or positive feedback with an 
enzymatic deactivation, such as is often found in signaling systems, 
the effects of a load are equally dramatic. We show that in a simple 
model of Ras activation, adding a small concentration of Raf 
molecules changes the bifurcation diagram of the signaling circuit 
and can completely abrogate the bistability in the system. While 
we have chosen a specific example of Ras activation, our simplified 
model, with an autocatalytic forward reaction and an enzymatic 
backward reaction is a minimal model for a many positive 
feedback switches. The change in the bifurcation diagram arises 
from the conservation laws that couple the concentration of the 
load with the concentrations of the proteins in the upstream module. 
Given the sensitivity of non-linear dynamical systems to initial 
conditions, it should probably be expected that many, if not all, 
positive feedback based switches that operate at the short timescales 
of signal transduction, and therefore must possess these conservation 
laws, should exhibit this sensitivity to the effect of a load. 

Our results throw up an interesting puzzle for quantitative 
biologists. In many natural signal transduction systems such as the 
MAPK cascade, the concentration of the output of a bistable 
switch is quite comparable to the concentration of the load, thus 
significant changes in load concentrations could have dramatic 
effects on the behavior of the switch. However it has also been 
shown that there is a significant cell to cell variability in protein 
concentrations [49]. How do cells ensure that positive feedback 
based switches such as the Ras switch continue to operate robustiy 
in the bistable regime? Additional regulatory mechanisms 
involving feedback between the load and its partner protein may 
exist that confer robustness to the qualitative behavior of the 
biochemical switch. Arguably some of the bells and whistles of 
natural protein networks that are often disregarded when 
analyzing the network may in fact be performing this role. In 
other words, self-assembled switches have to be complex! In this 



PLOS Computational Biology | www.ploscompbiol.org 



12 



IVIarch2014 | Volume 10 | Issue 3 | el 003533 



Loads Bias Biological Switches 



context it is worth mentioning that it has been persuasively argued 
[50,51] that some biological circuits maintain robustness of "fold- 
changi;' behavior rather than absolute levels of protein concen- 
tration. It is possible that additional protein-protein interactions 
that couple concentrations of loads with output proteins may end 
up in performing this function. Another significant factor that 
needs consideration is the role of spatial segregation in producing 
feedback from the downstream module to the upstream one. In 
fact it has been shown experimentally that MAPK substrates 
sequester activated MAPK in the nucleus, and thus protect it from 
cytoplasmic phosphatases. Changing the concentration of one- 
substrate therefore affects the concentration of activated MAPK 
[52]. 

Previous discussions of the effect of loads on the operation of 
circuits have suggested the use of insulators, that is circuit elements 
that insulate the upstream module from the downstream module 
[22]. The initial suggestions for building insulators in Ref [22] 
involved incorporating signal amplification along with negative 
feedback in the upstream circuit. Another way of insulating the 
circuit is to ensure that the demand of the load for its cognate 
repressor is never significant compared to the total amount of 
repressor. For a genetic switch therefore, a possible insulating 
mechanism is if the link to the downstream circuit is through a 
promoter. For example, consider making an AND gate from an 
output of the toggle switch. This can be done by inserting a 
constitutively produced protein Y that binds to Rl such that the 
complex is a transcription factor for another protein, say Z. Thus 
there is an AND relationship between the two inputs, Y and Rl 
and the output Z. To offset the effect of load induced modulation 
of the dynamics of Rl, an additional step can be inserted such that 
Rl first binds to the promoter region of another gene that codes 
for protein X and activates its transcription, and it is the protein X, 
rather than Rl, that can bind to Y and activate production of Z. 
The advantage of adding this extra step is that the concentration of 
the promoter for X is very small compared to the concentration of 
Rl, and therefore load induced modulation of the upstream toggle 
can be kept at a minimum. Note however that this cannot be done 
without the additional cost of the time delay required for the 
transcription and translation of X. 

As can be seen, any additional step or series of steps that can 
amplify a weak signal can act as an insulator. Another standard 
example of an amplifying circuit is a phosphorylation cascade 
which is especially relevant when considering Ras activation since 
it directly leads to the MAPK phosphorylation cascade. Phos- 
phorylation cascades are also very fast, and therefore do not face 
the additional time delays of an additional transcriptional step. 
From the point of view of .synthetic circuit design, the insulating 
mechanism here could be constructed by designing a weak binding 
affinity of Ras (or the synthetic protein that plays that role) for Raf 
(or the equivalent protein). The bound complex then catalyzes a 
phosphorylation cascade that ends by connecting to the down- 
stream circuit. 

Note that this method of insulation does not have the same time 
delay costs as the additional transcription steps. However it does 
come with the metabolic costs of having to produce large amounts 
of proteins that are essentially serving no useful physiological 
purpose for the cell. This cost could be relevant in some synthetic 
biology applications, and certainly needs to be evaluated during 
circuit design. It has been shown in the context of phosphorylation 
cycles that insulation always carries a metabolic cost, and in 
general better insulation carries a greater metabolic cost [53]. 

The existence of the MAPK phosphor)'lation cascade however 
begs the question whether it serves the purpose of insulation of the 
upstream Ras circuit from the downstream circuit. While it is not 



possible to answer this intriguing question without further 
experiments, it does appear that the Ras-Raf complex is present 
is quite large numbers on activated cells. This would suggest that 
insulation is not the function for which the cascade may have 
evolved. Our own analysis of the genetic toggle switch with the 
positive feedback motif suggests that Nature may prefer more 
complicated forms of regulation that balance the difiFerent 
components of the circuit. However there is no reason why both 
methods cannot be utilized. To our mind this is a very exciting 
question that requires more attention from experimentalists and 
theorists alike. 

It should also be noted that due to non-specific binding of 
transcription factors with DNA as well as between proteins, every 
circuit in the cell, real or synthetic, operates in the presence of a 
load. Variability in the functioning of circuits that are seen when 
transferring synthetic circuits between species, or even in different 
cells, may be a result of not only dififerences in basic protein 
concentrations, but also of this undervalued but nevertheless 
tangible load. Based on this reasoning we predict that some of the 
host-dependent effects that complicate synthetic biology, i.e. a 
synthetic circuit that works in one organism not performing well in 
another, are in fact due to changes in the intrinsic load due to non- 
specific binding when changing hosts. 

Our analysis underscores the importance of incorporating loads 
when simulating models of switches in natural and synthetic 
systems. Mathematical analysis of switch-like motifs therefore 
would do well to at least include a load on their output proteins, in 
order to incorporate the possible effects of load induced modu- 
lation on the circuit. 

Supporting Information 

Figure SI Surface plots showing response times of the 
simple genetic toggle switch with changes in load (L) and 

changes in the dissociation constant (Kd) of binding with 

load. The units of L and Kd are (m()lc('ulc's/|j.m' ). The z-axis 
measures the response time* indicated in the title. "Same Side 
Rise" and "Same Side Decay" refers to the rise time and decay 
time when the load is on the same side as the repressor whose 
concentration is increasing. "Opposite Side Rise" and "Opposite 
Side Decay" refers to the rise time and the decay time when the 
load is on the other side of the repressor whose concentration is 
increasing. "Both Sides Rise" or decay refer to the rise and decay 
times when a load is present on both sides (symmetrically). The 
plot shows that at every Kd, the relation between the response 
time and load is approximately linear. The response time is largest 
for the case of "Both Sides Rise" followed by "Opposite Side 
Rise". The response time is also non-monotonic with respect to 
the Kd for a given load, and is maximized at intermediate values 
of Kd. 
(TIF) 

Figure S2 Time plot of switching of the simple toggle 
switch with a load on Repressor 1, at three different 
values of the dissociation constant. In all three cases the 
system is switched by providing 150 molecules/nm^ of an inducer 
at 1000 minutes. The inducer stays constant at that value and is 
not shown in the plots. The left panel has a very high dissociation 
constant (Kd= 1000 molecules/(im'') of binding between the load 
and the repressor, due to which the load has a minimal effect on 
the system. The middle panel has an intermediate value (Kd= 1 
molecules/nm^) because of which the load acts as a dynamic sink 
by releasing Repressor 1 and slowing the switching. The right 
panel shows the effect of a small dissociation constant (Kd= 10 
molecules/nm^ At such strong binding affinities, all of the load is 
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always bound to Repressor 1 . Thus the load has minimal effect on 
the switching dynamics. In all cases total load concentration is 100 
molecules/nm^. 
(TIF) 

Figure S3 Effects of a dynamic load on dynamics of a 
symmetric toggle switch. (A). The time taken to reach 90% of 
maximum value for the protein undergoing a low-to-high 
transition as a function of the equilibrium constant of a dynamic 
load. Normalized time is a unit-less number defined by the 
transition time (rise or decay) of the system at a given loading 
condition divided by the transition time (rise or decay) of an 
unloaded system. (B). The time taken for the concentration of the 
protein undergoing a high-to-low transition to reach 10% of its 
maximum value. 
(TIF) 

Figure S4 Stochastic time trace and the probability 
distribution function of repressor concentrations for the 
large volume simulations. (A). Comparison of time traces of 
the stochastic simulations of the simple toggle switch with basal 
parameters (top panel) and a larger volume (bottom panel). The 
average molecule number is about 5 times greater, and the 
number of transitions are significantly fewer. (B). The probability 
distribution function of the genetic toggle switch with the larger 
molecular number without (left) and with (right) a load. The effect 
of a load on Rl is qualitatively the same for this system as for the 
smaller system. Since transitions are slower the data are more 
uneven for this simulation. 
(TIF) 

Figure S5 Transition times in a genetic toggle switch 
with a positive feedback moiety. In all cases the strength of 
the positive feedback (denoted here by P instead of p) is 3.5 on 
either Repressor 1 (Rl) or Repressor 2 (R2). Top Left: Rise tim(; - 
time to transition INTO state Rl with the positive feedback on 
Rl. Note that the rise time is larger at nonzero loads when the 
load is on R2 or when the load is on both sides, in agreement 
with the simple toggle switch. Top Right: Rise time - time to 
transition INTO state Rl with the positive feedback on R2. 
Bottom Left: Decay time - time to transition OUT OF state Rl 
with the positive feedback on Rl. Bottom Right: Decay time - 
time to transition OUT OF state Rl with the positive feedback 
on R2. 
(TIF) 

Figure S6 Probability distribution functions of repres- 
sor concentrations for the toggle with a positive 
feedback moiety. The left panel shows that when p = 0, the 
switch is balanced evenly. As p increases, the side of the switch 
with the positive feedback becomes more and more prominent, at 
the expense of the other side. When p = 5, the system spends most 
of its time in one state. 
(TIF) 

Figure S7 Bifurcation diagram of the Ras switch with 
different levels of Raf (load) on the system for the model 
with Pseudo Steady State Assumption (PSSA). The total 
number of SOS in the simulation box is used as the parameter 
being tuned, which \;u-ics from 0 to 1000. For Raf = 0, Raf= 10 
and Raf = 30, there are two bifurcations points as SOS is 
increased. In the first bifurcation a new high valued stable steady 
state appears along with the low valued stable steady state. In the 
second bifurcation, the low valued stable state disappears leaving 
behind only the high valued state. The dotted line marks the 
unstable steady state that also comes into existence in the bistable 
region. As total Raf increases, the two bifurcations approach each 



other. When Raf = 50, the system has lost both of its bifurcations 
and is characterized by a single stable steady state at all values of 
Raf 
(TIF) 

Figure S8 Bifurcation diagram of the Ras activation 
model based on Law of Mass Action (LMA) . Here the total 

number of Raf molecules (Rafi-) is the primary parameter being 
varied. Without Raf, the Ras activation system is bistable as 
reported. With increasing Rafx, the "high" stable steady state 
branch comes closer with the unstable steady state branch and 
both are eliminated after a threshold of Rafx. A monostable region 
is maintained beyond the threshold. 
(TIF) 

Figure S9 Bifurcation diagram of the Ras activation 
PSSA model with total number of Raf molecules (Raff) 
as the primary parameter. Without Raf, the Ras activation 
system is bistable as reported. With increasing Rafp, the "high" 
stable steady state branch comes closer with the unstable steady 
state branch and both are eliminated after a threshold of Rafx. A 
monostable region is maintained beyond the threshold. 
(TIF) 

Figure SIO Parameter sensitivity of the bistability of 

Ras switch to changes in k„£jg. Increase of k„fi^; results in 
leftward shifts of both stable fixed points, increase in the bistable 
regime and increase in maximal RasGTP activation level (Green 
Line) when compared to baseline with original value (Blue Line). 
Decrease of k^g^ (Red Line) results in right shift of both limit 
points, increase in unstable bistable regime and decrease in 
maximal RasGTP activation level. Qualitative features of bis- 
tabUity are maintained. 
(TIF) 

Figure Sll Parameter sensitivity of the bistability of 
Ras switch to changes in k„„s. Increase of kon6 (Green Line) 
results in right shift of both limit points, increase in unstable 
bistable regime and decrease in maximal RasGTP activation level 
when compared to baseline original value (Blue Line). Decrease of 
k<>n6 (Red Line) results in left shifts of both limit points, increase in 
bistable regime and increase in maximal RasGTP activation level. 
Qualitative features of bistability are maintained. 
(TIF) 

Figure S12 Comparison between bifurcation diagrams 
of toy genetic toggle switch with and without protection 
of repressor degradation when bound with promoters. If 

the protection is not included (Blue Line), a minor increase in the 
bistable region can be observed with right shift of upper limit point 
and left shift of louver limit point compared to the case with 
protection assumed (Red Line). Note that this is not the same as 
degradation after being bound with the load. 
(TIF) 

Figure S13 Bifurcation diagram of the genetic toggle 
switch with positive feedback loop on one side after 
removal of the protection assumption. The left panel shows 
the bifurcation diagram when the load is added symmetrically to 

both sides. Without load molecule, the toggle switch is bistable as 
predicted. With the increase in Lx, the unstable steady state and 
the "low" stable steady state come closer and meet at certain 
threshold. The value of "high" stable steady state decreases with 
increase in Lx. Beyond the threshold, the toggle switch becomes 
monostable. The right panel shows the effect of just adding a load 
to Rl. In this case the high state of Rl approaches the unstable 
steady state, and annihilates itself The system jumps to the low 
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stable state, which is equivalent to the "high" state of the other 

repressor. 

(TIF) 

Figure S14 Bifurcation diagram of the Ras activation 
model when Ras can degrade when bound with Raf. As 

the number of Raf molecules increase, the bistable region 
decreases. However unlike the case with no protection, the curve 
moves to the left. When Raf molecules increase by a large amount, 
bistabUity is abrogated. 
(TIF) 

Figure S15 Transition times for various k'o„ and k'„£F 
values plotted as a function of load for the basic toggle 
switch. Even if the binding-unbinding rates are slower or much 
faster than protein decay rates, the load-transition time relation- 
ship stays linear. A, C, E, G, I, K, M, O, Q^and S .show the rise 
time. B, D, F, H, J, L, N, P, R, and T show decay time. (A,B) 
k'„n = 4, k'„ff=0.5, Kd = 0.12."5. (C,D) k'„„=l()', k'„ff=0.5, 
Kd = 0.05. (E,F) k'o„ = 4, k'off=4, Kd=l. (G,H) k'o„=10, 
kV=10, Kd=l. (IJ) k'„„ = 4, kV=20, Kd = 5. (K,L) 
k'„„=10, k'„ff=,50, Kd = 5. (M,N) k'„„ = 4, k'„ff=40, Kd=10. 
(0,P) k'on= 10, k'„ff= 100, Kd= 10. (Q^R) k'on = 40, k'off=400, 
Kd= 10. (S,T) k'o„= 100, k'off = 1000, Kd= 10. 
(TIF) 

Figure SI 6 Probability distributions of repressor con- 
centrations for various values of k '„„ and k'„fr for the 
basic toggle switch. Even when the binding-unbinding with the 
load is several times faster than protein decay rates, the basic 

phenomena discussed in the paper remains unchanged. (A) 
k'on = 50, k'off=500 (B) k'o„ = 500 k'„ff= 500 (C) k'on = 500 
k'off=5000. 
(TIF) 

Figure S17 Bistability of the toggle switch with positive 

feedback. A bifurcation diagram of the simple toggle switch 
with a positive feedback moiety on one side, with respect to the 
parameter p that measures the strength of the positive feedback. 
Only the concentration of Rl is shown for simplicity. The 
switch remains bistable till p becomes larger than a little over 
200. 
(TIF) 

Table SI Slopes of linear fits to rise and decay time 
with various values of Kog-, K„„ and p. The first column 
reports the values of the dissociation constant (Kd = Koff/Kon) 
and the kinetic constants of the binding of Repressor 1 , 2 or the 
value for P, which represents promoter strength. The other 
columns report the slopes of the linear fits of the various rise 
times and dcciu" times. In all ca.ses the fits have- high R-s(|uared 
values (>0.95). Intercept is 1, as the slopes are normalized to 
the un-loaded transition time. For K^j we change the para- 
meters by two orders of magnitude in both directions to show 
that the linear relation is robust despite these changes. Note 
that the relation between rise time or decay time and the 
binding constant is non-monotonic. Units are as reported in the 
text. 
(DOC) 

Table S2 Exponential Fits of the amount of inducer 
required to transition states as a function of load. The 

basic genetic toggle switch switch was toggled to its other state by 
production of the other repressor protein by an inducer, given 
here as a bolus with a decay rate as shown. The size of the bolus 
was increased until the state changed. This was repeated at 
different levels of load and the minimum size of the bolus 



required was fit by an exponential function of the load. The fits 
are shown here, along with their R-squared values. "Load 
applied to the opposite side" means switching from a state 
without a load to a state with a load. "Load applied to the same 
side" means switching from a state with a load to a state without 
a load. 
(DOC) 

Table S3 Exponential Fits of the amount of inducer 
required to transition states as a function of load, in the 
case of induction by repression. The switch was toggled to its 
other state by repression of the current state by an external 

molecule, given to the system as a bolus with a decay rate as 
shown. The size of the bolus was increased until the state changed. 
This was repeated at different levels of load and the minimum size 
of the bolus required was fit by an exponential function of the load. 
The fits are shown here, along with their R-squared value. Thus 
the inducer required depends exponentially on the load in both the 
methods of induction. "Load applied to the opposite side" means 
switching from a state without a load to a state with a load. "Load 
applied to the same side" means switching from a state with a load 
to a state without a load. 
(DOC) 

Table S4 Slopes of linear fits to rise and decay time 
with a dynamic load, with varying values of load decay 
rate K^, load binding rates K^^ and and constant 

Ki. The first four columns report the values of the various 
parameters. The other columns report the slopes of the linear fits 
of the various rise times and decay times. In most cases the fits 
have high R-squared values (>0.9.'i). The two exceptions are 
>0.90 and starred. Intercept is 1, as the slopes are normalized to 
the un-loaded transition time. Note that for all cases, the 
relationship between load (expressed here as Keq = Ki,/Kd) and 
transition time is a positive linear relationship. 
(DOC) 

Table S5 Rate expressions used for the stochastic 
simulations of the genetic toggle switch. The rate 

expressions used for the stochastic simulation of the toggle switch 

along with the description of the reaction are listed. 

(DOC) 

Table S6 List of reactions in the minimal model of Ras 
activation. The reactions in the minimal model of Ras 
activation, jilong with the labels of the corresponding rate 
constants are shown. Parameters used in the simulations are given 
in Table S7. 
(DOCX) 

Table S7 Kinetic rate parameters used for the simula- 
tions of the Ras model. Here the numbers in the subscript of 

the rate constants in the "Constant" column refer to the reactions 
shown in the corresponding row of Supplementary Table S6. The 
meaning of the rate constants are as follows: ^^an refers to the on- 
rate, k^gk the off rate and k^at is the catalytic rate. The sources for 
the rates are as shown in the last column. 
(DOC) 

Table S8 List of reactions in the toy model of genetic 
toggle switch. The reactions in the toy model of the genetic 

toggle switch, discussed in Supplementary Text SI section 3.1 are 
listed. The description of the various chemical species in the 
reactions are also provided in the Supplementary Text SI. 
(DOCX) 

Text SI Supporting Information including derivation 
and analysis of toggle switch and Ras models, param- 
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eters used for the simulation and their sources, results 
of the parameter sensitivity analysis, details of the effect 
of a dynamic load on the genetic toggle switch, results 
for the Ras model with Michealis-Menton kinetics, 
results for the models without the protection assump- 
tion. 
(PDF) 
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